int k=2, lb=0;
#include "listPlex.h"
namespace ListPlex
{
    template<typename VType>
    struct PlexEmitor{
        graph<VType> g;
        long long cnt=0;
        // cilk::reducer_opadd<long long> cnt;
        std::function<VType(VType)> sub2ori;
        template<typename Fn>
        PlexEmitor(graph<VType> &_g,Fn&& f):g(_g), sub2ori(std::forward<Fn>(f)), cnt(0){}

        void emitPlex(VType *list, int len){
            //TODO: atomic or thread-own
            // utils::CAS(&cnt,cnt,cnt+1);
            cnt++;
        }
    };

    template <typename VType = int>
    struct DecomposeMaxChecker
    {
        graph<VType> g;//global graph
        int nSub;
        VType *sub2ori; 
        VType *ori2sub;
        VType vstart;
        VType *leadList;
        int idx;
        std::function<bool(VType, VType)> predFunc;
        
        DecomposeMaxChecker(graph<VType> &_g,int _nsub, VType *_sub2original, VType *_original2sub, VType _start, VType *_preds, int _len, std::function<bool(VType, VType)> &&_p) 
            : g(_g),nSub(_nsub), sub2ori(_sub2original), ori2sub(_original2sub), vstart(_start), leadList(_preds), idx(_len), predFunc(_p){}

        ~DecomposeMaxChecker(){
            delete[] ori2sub;
            delete[] sub2ori; 
        }

        //Counter the number of adjacent vertices of oriU in the growing plex      
        bool canGlobalAdd(VType oriU, int szplex, MBitSet32 &plex){
            //para_for is not good
            int cnt = 0;
            for (int j = 0; j < g.V[oriU].degree; j++) {
                VType nei = g.V[oriU].Neighbors[j];
                VType neiSub = ori2sub[nei];
                //vertex neiSub is in the kplex
                if (neiSub < nSub &&  plex.test(neiSub)) {
                    cnt++;
                }
            }
            if (cnt + k >= szplex + 1)
                return 1; //can add
            return 0;	//can not add
        }

        
        //Check if the vertices in list in a globaly maximal k-plex
        //i.e., maximality : all vertices in Preds cannot form a larger k-plex with list.
        bool isMaximal(VType *list, int szplex, int* neiInP)
        {   
            // saturated vertices;
            assert(szplex < PSIZE_LIMIT);
            VType sat[PSIZE_LIMIT]; //saturated vertices in plex
            int nsat = 0;
            MBitSet32 plex(nSub); //mark a plex vertex in sub g;
            bool max = true;

            for (int i = 0; i < szplex; i++){
                int v = list[i];
                plex.set(v);
                if (neiInP[v] + k == szplex){
                    sat[nsat++] = v;
                }
            }            
            if (nsat > 0){
                VType vfirst = sub2ori[sat[0]];
                VType *commons = new VType[g.V[vfirst].degree];                             
                copy(g.V[vfirst].Neighbors, g.V[vfirst].Neighbors + g.V[vfirst].degree, commons);
                int szcom = g.V[vfirst].degree;
                for (int i = 1; i < nsat; i++){
                    VType u = sub2ori[sat[i]];
                    szcom = utils::interSection(g.V[u].Neighbors, g.V[u].degree, commons, szcom, commons);                    
                }
                //(commons,szcom) contains all common adajcent vertices of saturated vertices                
                for (int i = 0; i < szcom; i++){
                    int u = commons[i];
                    if (predFunc(u, vstart) && canGlobalAdd(u, szplex, plex)){ // not maximal because we can add u to the solution
                        max = false;
                        break;
                    }
                }   
                delete[] commons;             
            }else{//no saturated vertex
                for (int i = 0; i < idx; i++){
                    VType upred = leadList[i];
                    if (canGlobalAdd(upred, szplex, plex)){
                        max = false;
                        break;
                    }
                }
            }
            //return true;
            return max;
        }           
    };

    template<typename Pred>
    _seq<VType> decompose(Graph &g,int vstart, Pred &&pred){
        int *dist = new int[g.n]; // the dist to vstart 
        memset(dist, 0, sizeof(int) * g.n);
        //parallel_for(int j = 0; j < g.n; j++) dist[j] = 0;

        vector<VType> neibors_1hop;
        vector<VType> neibors_1a2hop;

        //collect 1 hop neigbors
        dist[vstart]  = 1;
        for(int j = 0; j < g.V[vstart].degree; j++){
            VType nei = g.V[vstart].Neighbors[j];                    
            if (!dist[nei] && pred(vstart, nei)){// vstart pred nei
                dist[nei] = 1;
                neibors_1hop.emplace_back(nei);
                neibors_1a2hop.emplace_back(nei);
            }
        }
        
        //collect 2-hop neigbors
        for (auto nei1: neibors_1hop){
            for (int j = 0; j < g.V[nei1].degree; j++){
                VType nei2 = g.V[nei1].Neighbors[j];
                if (!dist[nei2] && pred(vstart, nei2)){
                    dist[nei2] = 2;
                    neibors_1a2hop.emplace_back(nei2);
                }
            }
        }        
        prompt("[%d] 2-hop neighbors %d\n", vstart, (int)neibors_1a2hop.size());
       
        //pruning the  vertices
        //the common neigbors between v and a 2-hop neighbor must be above certain threshold
        // int nRed = 0;
        cilk::reducer_opadd<int> nRed(0);
        int sz_neibors_1a2hop = neibors_1a2hop.size();
        parallel_for (int j = 0; j < sz_neibors_1a2hop;j++)
        {
            int nei = neibors_1a2hop[j];
            int common = 0;
            int j1=0, j2=0;
            while (j1 != g.V[vstart].degree && j2 != g.V[nei].degree){
                //NOTE: vertex in Neighbors are sorted in ascending order
                if (g.V[vstart].Neighbors[j1] < g.V[nei].Neighbors[j2]){
                    j1++;
                }else if (g.V[vstart].Neighbors[j1] > g.V[nei].Neighbors[j2]){
                    j2++;
                }else{
                    //A common neighbor.
                    if (dist[g.V[vstart].Neighbors[j1]]){ 
                        common++;
                    }
                    j1++,j2++;
                }
            }
            if (dist[nei] == 1){//1-hop neibor
                if (common + 2 * k < lb) {
                    dist[nei] = 0; //remove nei
                    nRed++;
                }
            }else{ //dist = 2
                if (common + 2 * k - 2 < lb) {
                    dist[nei] = 0; //remove nei
                    nRed++;
                }
            }
        }
        int nRedV = nRed.get_value();
        // int nRedV = nRed;
        prompt("Reduce %d, remaining %d\n", nRedV, (int)neibors_1a2hop.size() - nRedV);
        
        //VType *block = sequence::pack(neibors_1a2hop.data(), dist, g.n);
        VType *block = newA(VType, neibors_1a2hop.size() - nRedV + 1);        
        int szblk = 0;
        block[szblk++] = vstart;
        for (auto v: neibors_1a2hop){
            if (dist[v] != 0){
                assert(v != vstart);
                block[szblk++] = v;
            }
        }
            
        assert(szblk == neibors_1a2hop.size() - nRedV + 1)  ;
        delete[] dist; 
        return _seq<VType>(block, szblk);
       //return _seq<VType>(nullptr, 0);
    }

    Graph buildSubgraph(const Graph &g,_seq<VType> block, VType* map_org2sub, VType *map_sub2org){            
        int nv = block.n;
        VType *deg = new VType[nv];
        memset(deg, 0, sizeof(VType) * nv);
        VType ne = 0;
        //int *offset = new int[nv];        
        for (int i = 0; i < block.n; i++){
            VType original = block.A[i];
            //int deg = 0;            
            //printf("[%d]:",original);
            for (int j = 0; j < g.V[original].degree; j++){
                VType nei = g.V[original].Neighbors[j];
                if (map_org2sub[nei] < nv){ // in block
                    //printf("%d ",nei);
                    deg[i]++;
                    ne++;
                }
            }
            //printf("\n");
        }
        //printf("ne:%d\n", ne);
        //build the graph
        VType *in = newA(VType, nv + ne + 2);
        in[0] = nv;
        in[1] = ne;
        VType* offset = in + 2;
        VType* edges = in + 2 + nv;
        //use SCAN inclusive
        for (int i = 0; i < nv; i++){
            offset[i] = (i==0)? 0 : offset[i-1] + deg[i-1];           
        }               
        vertex<VType> *vertices = newA(vertex<VType>, nv);
        parallel_for(int i = 0; i < nv; i++){
            vertices[i].degree = deg[i];
            vertices[i].Neighbors = edges + offset[i];
            //TODO add edges;            
            VType original = block.A[i];
            int cnt = 0;            
            for (int j = 0; j < nv; j++){
                VType nei = map_sub2org[j];
                bool isnei = binary_search(g.V[original].Neighbors, g.V[original].Neighbors + g.V[original].degree, nei);
                if (isnei){
                    vertices[i].Neighbors[cnt++] = j;
                }
            }
            assert(cnt == deg[i]);            
        }        
        delete[] deg;
        return Graph(vertices, nv, ne, in);   
    }

    void decomposableSearch(Graph &g){
        VType *dseq = new VType[g.n];
        VType *core = new VType[g.n];
        VType *dpos = new VType[g.n];
        bool *mark = new bool[g.n]; 
        auto pred = [dpos](VType v1 ,VType v2){
            return dpos[v1] < dpos[v2];
        };

        Timer watch;
        watch.start();

        VType maxcore = degeneracyOrder<VType>(g, dseq, core, dpos);        
        prompt("Degeneracy. Max core %d. List: ", maxcore);
        
        //Collect all the possible leading vertex, i.e. core[v] + k >= lb. store in leadList
        #pragma cilk grainsize = 128
        parallel_for (int i = 0; i < g.n; i++ ){
            mark[i] = core[dseq[i]] + k >= lb? true: false;            
        }
        _seq<VType> leadList = sequence::pack(dseq, mark, g.n);            
        prompt("\nNumber of possible vertices %ld\n. Now decompose. \n", leadList.n);

        cilk::reducer_opadd<intT> validblk(0);
        cilk::reducer_opadd<long> cnt_maxplex(0);
        // intT validblk=0;
        // long cnt_maxplex=0;

        parallel_for(int idx = 0; idx < leadList.n; idx++){ //parallelizable
            // Search with leading vertex vstart;
            VType vstart = leadList.A[idx];
            prompt("Leading vertex %d\n",vstart);

            // decompose the graph with leading vertex vstart
            _seq<VType> blk = decompose(g, vstart, pred);
            
            if (blk.n < lb) {
                blk.del();            
            }else{            
                //build subgraph 
                int *map_org2sub = new int[g.n]; //map from orinial vertex to subgraph vertex
                int *map_sub2org = new int[blk.n];                
                for (int j = 0; j < g.n; j++) map_org2sub[j] = g.n;
                //memset(map_org2sub, 0x3f, sizeof(int) * g.n); // initialize
                parallel_for (int isub = 0; isub < blk.n; isub++){
                    VType original = blk.A[isub];
                    map_sub2org[isub] = original;
                    map_org2sub[original] =  isub;
                }
                //build graph
                Graph subg = buildSubgraph(g, blk, map_org2sub, map_sub2org);
                //prompt("subgraph %d: n:%d m:%d. Vertex:", validblk, subg.n, subg.m);               
                for (int j = 0; j < blk.n; j++) 
                    shortprompt("%d ", blk.A[j]);
                shortprompt("\n");                

                ///build search context
                int *neiInP = new int[subg.n];
                int *neiInG = new int[subg.n];
                PlexStack *plex = new PlexStack(PSIZE_LIMIT);
                VtxSet *excl = new VtxSet(subg.n);
                VtxSet *cand = new VtxSet(subg.n);
                for (int j = 0; j < subg.n; j++)
                {
                    neiInP[j] = 0;
                    neiInG[j] = subg.V[j].degree;
                    cand->add(j);
                }
                //graph<VType> &_g, int _k, int _nsub, VType *_sub2original, VType *_original2sub, VType _start, VType *_preds, int _len, Pred &_p                
                DecomposeMaxChecker<VType> maxchecker(g,subg.n, map_sub2org, map_org2sub, vstart, leadList.A, idx, pred);
                auto vtxmapfunc = [&] (VType v){return map_sub2org[v];};
                PlexEmitor<VType> emitor(g,vtxmapfunc);
                KplexListor<DecomposeMaxChecker<VType>, PlexEmitor<VType>,VType> listor(
                        subg,  &maxchecker, &emitor, plex, cand, excl, neiInP, neiInG);
                listor.candToPlex(0)->listByCase();//start from the first vertex
                // listor.branchBySelection(0);//start from the first vertex

                // cnt_maxplex += emitor.cnt.get_value();
                cnt_maxplex += emitor.cnt;
                validblk++;           
                blk.del();//release blk  
                subg.del();
            }
        }
        double durat = watch.stop();
        printf("Total valid block %d\n", validblk.get_value());
        // printf("Total valid block %d\n", validblk);
        // printf("Number of %d-plex: %ld\n", k, cnt_maxplex);
        printf("Number of %d-plex: %ld\n", k, cnt_maxplex.get_value());
        printf("Running time %.3f\n",durat);

        delete[] dseq;
        delete[] core;
        delete[] dpos;
        delete[] mark;
        leadList.del();
    }
}
int main(int argc, char **argv)
{
    graph<intT> g = ListPlex::readBinaryGraph(argv[1]); 
    if (argc == 3){       
        k = atoi(argv[2]);                        
        ListPlex::decomposableSearch(g);
    }else if(argc == 4){
        k = atoi(argv[2]);
        lb = atoi(argv[3]);
        ListPlex::decomposableSearch(g);
    }else {
        fprintf(stderr, "usage: enplex <filename> <k> [lb=0] [cut-off time=100(sec)]\n");
    }
    return 0;
}